MODELACIÓN DEL AGUA SUBTERRÁNEA A ESCALA 
REGIONAL CON REFINAMIENTO LOCAL DE LA 


MALLA. PLANTEAMIENTO Y VALIDACIÓN 
DEL ALGORITMO 


e Eric Cabrera-Estupiñán e 
e Armando Hernández-Valdés e 
Instituto Superior Politécnico José Antonio Echeverría, Cuba 


Resumen 


La simulación del flujo del agua subterránea a escala regional no permite 
reproducir los problemas locales que provocan descensos significativos de los 
niveles del agua en las obras de captación. A pesar de tal afirmación, este tipo 
de modelos se ha utilizado para evaluar el comportamiento hidrodinámico 
del sistema acuífero y obras de captación a esa escala. En el presente trabajo se 
realiza la validación de un algoritmo que permite, a partir de un modelo regional, 
reproducir el comportamiento local de los niveles próximos a pozos de bombeo 
que por su importancia requieran de un control sistemático para ser objeto de 
operación en tiempo real mediante técnicas de control automático. A través de la 
validación con soluciones analíticas, se demuestra la importancia del refinamiento 
de la discretización, incluso empleando el triángulo cuadrático en el Método de 
los Elementos Finitos, que es el utilizado en la tecnología AQUIMPE, ahora en 
su versión sobre Windows. El algoritmo propuesto permite reproducir los efectos 
de la no linealidad del flujo en las proximidades de los pozos de bombeo y los 
efectos de la penetración parcial de estos, mediante modificaciones en los valores 
de transmisividad de los elementos cercanos al pozo de bombeo de acuerdo con 
los rangos de distancia utilizados en el algoritmo numérico. 


Palabras clave: modelación matemática, acuíferos, recursos hidráulicos, sistemas 
de información geográfica, método de elemento finito, refinamiento de malla. 


Introducción de cómputo; por otra parte, se pierde la 
posibilidad de utilizar elementos volumétricos 


En la modelación de los sistemas acuíferos representativos de medios porosos equiva- 


a escala regional se establecen criterios para 
definir el tamaño de los elementos de acuerdo 
con el método numérico empleado y los 
objetivos del modelo (Reilly y Harbaugh, 2004; 
Hernández et al., 2001; Jagelke y Barthel, 2005). 
Todos estos autores coinciden en que con mallas 
gruesas no se puede pretender representar los 
efectos locales de pozos y campo de pozos; pero, 
por otra parte, mallas muy finas complican el 
proceso de calibración y al no disponerse en 
muchas ocasiones de bases de datos que la 
justifiquen, se incrementa significativamente 
el costo de la simulación en cuanto a tiempo 


lentes en zonas de fisuras y fracturas, como son 
los medios cársicos. 

Hasta la fecha se han utilizado los modelos 
de simulación regional para evaluar políticas 
existen 
cargas 
simuladas y las correspondientes originadas 


de explotación, conociendo que 


diferencias significativas entre las 
en las obras de captación. La utilización del 
Método de los Elementos Finitos (MEF), con el 
triángulo cuadrático empleado en la tecnología 
AQUIMPE, permite reducir el número de 
elementos, ya que dentro de cada triángulo 


se ajusta la superficie solución, pudiéndose 
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determinar la carga piezométrica en cada 
uno de sus puntos mediante una superficie 
cuadrática aproximada, de forma tal que exista 
coincidencia en los seis nodos del triángulo 
(Martínez, 1989). Los elementos pequeños 
se emplean con este método en las zonas 
de mayor complejidad geológica, donde las 
propiedades hidráulicas y geométricas del 
acuífero tienen gran variación, lo cual se refleja 
en el comportamiento de las hidroisohipsas. 
También se utilizan elementos pequeños en 
zonas de intensa explotación o recarga, y 
donde se quiere aumentar la información sobre 
las respuestas del acuífero. En ambos casos, el 
tamaño del elemento no debe ser inferior a la 
base informativa disponible (Hernández et al., 
2001). 

Para considerar los efectos locales de los 
pozos de bombeo en la modelación regional 
se han utilizado diversos procedimientos. 
Según Anderson y Woessner (1992), al 
utilizar el Método de las Diferencias Finitas 
(MDP), Prickett (1967) y Trescott et al. (1976) 
propusieron emplear ecuaciones analíticas 
en función del tamaño y la forma de la celda 
para calcular las cargas en los pozos de 
bombeo. Si bien es cierto que con el Método 
de los Elementos Finitos se obtienen mejores 
aproximaciones de la carga cuando los nodos 
coinciden con pozos de bombeo, los autores 
discrepan de lo planteado por Anderson y 
Woessner (1992) en relación con que “no es 
necesario emplear fórmulas de corrección”, 
ya que siempre estará presente un Error de 
Refinamiento de la Discretización (ERD), 
como se demuestra en el presente trabajo. 

Otro procedimiento empleado para con- 
siderar los efectos locales consiste en un 
refinamiento de la discretización del modelo 
regional como la realizada con la tecnología 
AQUIMPE en el modelo CIRO, con la creación 
del modelo RUSPOLlI, del campo de pozos del 
acueducto que abastece a la ciudad de Ciego de 
Ávila, en Cuba (Hernández, 1991). En este caso, 
la simulación se realizó suministrándole las 
cargas del modelo regional a los contornos del 
modelo local, pero en ningún caso se pudieron 


simular los abatimientos que se producían en 
los pozos de bombeo. 

Más recientemente se han realizado aportes 
en el campo del refinamiento de mallas para el 
MDE, entre los que se pueden citar las técnicas 
de Refinamiento de Malla Gradual (GMR), 
Refinamiento de Malla Local (LGR) y Refina- 
miento de Malla Telescópica. Todas estas 
técnicas han sido utilizadas con el simulador 
MODFLOW (Leake y Claar, 1999) y otros 
modelos como ZOOMQ3D (Matthew et al., 
2006). En este campo se continúa trabajando 
intensamente y se brindan nuevos métodos 
para mejorar los acoples necesarios entre las 
mallas regionales (cuadrículas más gruesas) y 
las locales (cuadrículas más finas) (Scott et al., 
2006). 

En relación con el MEF, se pueden 
destacar los generadores de mallas flexibles 
que tiene el modelo FEFLOW, en donde la 
discretización en elementos finitos permite 
usar mallados complejos poco estructurados, 
con los cuales se logra una gran coincidencia 
en la representación de estructuras naturales 
mientras se manejan requerimientos como el 
tamaño del elemento, etcétera. Para grandes 
áreas de modelación, la generación de elemen- 
tos está soportada por sofisticados algoritmos 
de creación automática de la malla para 
garantizar un trabajo eficiente. Las mallas 
generadas automáticamente también tienen 
que ser adaptadas a las estructuras geográficas 
internas y hasta cierto punto locales, como son 
los ríos y pozos de extracción (DHI-WASY, 
2009). 

La nueva versión de AQUIMPE sobre 
Windows permite de manera directa simular las 
cargas que se generan en los alrededores de 
los pozos de bombeo sin desligar la influencia 
del modelo regional. Las particularidades 
de cada pozo en cuanto a su estructura, 
penetración parcial, pérdidas de carga por 
no linealidad del flujo, etcétera, pudieran ser 
calibradas modificando las propiedades de los 
elementos alrededor de cada pozo de bombeo 
de acuerdo con las cargas observadas en estos 
y en pozos satélites próximos a los mismos. 
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Estos modelos que integran los problemas 
locales al modelo regional permiten utilizar 
los datos de las propiedades obtenidas de la 
calibración regional y las asignaciones de la 
recarga y explotación difusa que influyen en 
el comportamiento de las obras de captación. 
La interferencia entre los pozos, fronteras 
hidrogeológicas y los efectos de las variaciones 
frecuentes en los caudales de explotación 
pudieran ser simulados y vinculados con 
la adquisición de datos en tiempo real para 
realizar modificaciones en las políticas de 
operación a corto plazo de los pozos y campo 
de pozos (Gómez, 2009). 

Estas ideas pueden materializarse en el 
presente trabajo también gracias a la creación 
de la herramienta AQTRIGEO que es una 
plantilla sobre la base de Auto CAD Map 3D 
2005, la cual tiene incluida bloques, capas, 
topologías, consultas a las topologías y macros 
en visual Basic for Application. Dicha plantilla 
es una herramienta SIG que sirve como 
plataforma para la producción y gestión de 
toda la base informativa de carácter espacial 
y de atributos que necesita AQUIMPE. 
Especialmente, este sistema permite la 
construcción por parte del especialista de la 


malla numérica para el modelo regional, la 
cual se realiza de forma semiautomatizada, ya 
que el criterio del modelador en la ubicación 
de los nodos y triángulos es importante y se 
le deja a éste, quedando solamente el proceso 
de numeración de los nodos y elementos de 
la malla de forma automática (ver Cabrera, 
2007). Además se realizó una nueva macro en 
AQTRIGEO, que permite captar la ubicación 
espacial de todos los pozos de explotación 
y realizarles un mallado fino a su alrededor 
de forma automática para posteriormente 
unir estas nuevas zonas entre ellas y con los 
triángulos del modelo regional (ver Cabrera y 
Escartín, 2008). 

En el presente trabajo se destaca el algo- 
ritmo de integración entre el modelo regional 
y el local, y su validación con las soluciones 
analíticas o teóricas. 


Modelo regional y modelos regional con local 
integrados 


En la figura 1 se muestra la discretización del 
modelo regional hipotético MODER 1, con 24 
elementos, representando un acuífero de 2 000 
m de ancho y 5 000 m de longitud, con 100 


5 000 m 
A 


Figura 1. Discretización del modelo regional hipotético MODER 1. 
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m de espesor, 50 m/día de conductividad 
hidráulica de Darcy K,, y coeficiente de 
almacenamiento E = 0.2. El nivel estático 
inicial es horizontal y con cargas conocidas 
fijas en sus fronteras izquierda y derecha, e 
impermeables en la superior e inferior. Un 
pozo bombea a caudal constante en un punto 
situado a 1 000 m de la frontera izquierda e 
igual distancia de las fronteras superior e 
inferior, justo en el centroide de la cuadrícula 
formada por los nodos [1,43, 45 y 3]. En este 
modelo, el pozo se concibe como un nodo 
principal de la malla. 

Para considerar los efectos locales de los 
pozos de bombeo en la modelación regional, 
se ha propuesto la generación de una malla 
alrededor de cada pozo de bombeo de interés, 
que parte desde el radio del pozo de bombeo 


con una progresión en distancias de los nodos 
principales medidos desde el centro de cada 
pozo y siguiendo una ley exponencial del tipo 
3", donde n toma los valores (-1, 0, 1 ,2 y 3). 
Finalmente, los nodos extremos se enlazan a 
nodos equivalentes de otros pozos o al modelo 
regional correspondiente, como se puede 
observar en la figura 2, con la discretización 
utilizada en MODEL 1. Aquí puede observarse 
que a diferencia del modelo MODER1, donde 
el pozo es simulado por un nodo del modelo, 
en este caso el pozo de bombeo está simulado 
por dos triángulos (34 y 35), que tienen 
aproximadamente la misma dimensión que el 
pozo real. 

Para evaluar el error por efecto del 
refinamiento de la discretización, se utilizó el 
modelo MODEL 2, que se muestra en la figura 


2 000 m 


5000 m 


126 129 


/ 


Figura 2. Discretización del modelo MODEL 1, con un conjunto de acercamientos graduales hasta el pozo. 
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Figura 3. Discretización del modelo MODEL 2, con un conjunto de acercamientos graduales hasta el pozo. 


3, de forma tal que una misma área de nueve 
hectáreas estaría discretizada en MODER 1 
con ocho elementos, MODEL 1 con 46 elemen- 
tos y MODEL 2 con 82 elementos. 


Simulación del flujo no lineal 


Entre los efectos locales que provocan mayores 
abatimientos en los pozos de bombeo que 
los calculados con régimen lineal de Darcy 
pueden estar los efectos de la presencia de un 
flujo no lineal y la penetración parcial; ambos 
fenómenos pueden ser simulados por una 


reducción de la transmisividad de Darcy en 
los elementos, según se aproximan al pozo de 
bombeo. 

Para este análisis se considera importante 
introducir las expresiones que determinan el 
abatimiento en cualquier zona del acuífero, 
producto del bombeo en un pozo: 


S=Sp+8, =0.183 2 los (7) 
Tb rE 


2 (1) 
Ñ : ] : 
217) r 
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2 
S =0.366 [e Q ] 1 (ta) 
Tp r 2nT7 ) r 


2.246T,+ 


1 
% 
De donde 7,= E ] y se conoce como 


radio de influencia. 

En las expresiones anteriores, el primer 
sumando de la derecha representa la 
componente darciana del abatimiento y el 
segundo término la componente turbulenta, 
siendo Q el caudal de bombeo; S, el abatimiento; 
T, = mK,, la transmisividad marciana; T, = 
mXK,, la transmisividad turbulenta, donde m 
es el espesor del acuífero; E, el coeficiente de 
almacenamiento; t, el tiempo de bombeo, r, 
la distancia desde el pozo de observación al 
centro del pozo de bombeo. 


De las expresiones anteriores y a partir 


de las conductividades hidráulicas darciana 


T Ea 
7 y turbulenta Ko= =p y por las 
siguientes ecuaciones: KE y Ki= Ea. 


donde v es la viscosidad cinemática del fluido 


y g la aceleración de la gravedad, se pudieran 
obtener las propiedades hidrogeológicas 
características de los acuíferos: permeabilidad 
intrínseca O geométrica (k) y rugosidad 
equivalente del medio poroso o fisurado (C). 

La diferencia de abatimientos entre dos 
puntos situados a distancias r, y r, del pozo 
de bombeo puede ser determinada a partir de 
la expresión general del abatimiento (Pérez, 
2001): 


S,-S,= AS,,= 0.366 o ls|7) 


D f 
E 
Q Y Y 
2 1 


La expresión anterior puede ser expresada 


Q) 


en función del radio de Darcy, definido por 
Pérez (2001), por la ecuación: 


A EA 
9017 7) 5 


Por lo que la ecuación (2) se transforma en: 


5,=5,= AS ¿0366 2 log | 2 
Tb 11 (4) 


Qt t=ñ 
Tp 40T 11) 


La diferencia de abatimientos anterior 
puede ser simulada por un medio heterogéneo 
equivalente de transmisividad T',, variable para 
cada rango de distancias 1, y 7,. 


S,-S,=AS,, =0.366 - log (=) (5) 


D 


Igualando las dos expresiones anteriores 
se obtiene la relación entre las dos transmisivi- 
dades darcianas, expresada por: 


La ecuación anterior permite obtener 
numéricamente por rangos de distancias las 
transmisividades de Darcy equivalentes, que 
producirían los mismos abatimientos que si 
se considerara un régimen no lineal de flujo 
a partir de la ecuación diferencial en régimen 
no lineal (Hernández, 1982), pero como medio 
homogéneo. 

En el cuadro 1 se presentan los valores 
del inverso del coeficiente a. hallados para 
un conjunto de relaciones (1, y r,) referidas al 
modelo MODEL 1 y para distintos radios de 
Darcy. 

De acuerdo con el rango de distancias 
anteriores, se asigna a los elementos que le 
correspondan a ellas una reducción porcentual 
de la conductividad hidráulica regional, 
según el valor estimado del radio de Darcy 
obtenido de las pruebas de bombeo o por un 
proceso de calibración. 
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Cuadro 1. Factores para afectar la transmisividad de Darcy a los elementos del modelo MODEL 1, 
según rango de distancias y radios de Darcy. 


Valores de (1/0) para distintos radios de Darcy 
r, (m) r, (m) MODEL 1 (0) 
r,=25 m r,=50 m r,= 80 m r,= 120 m 
0.25 il 1 + 0.108 r,, 0.27 0.156 0.104 0.072 
1 3 1 + 0.0304 r,, 0.568 0.397 0.291 0.215 
3 9 1 + OOO 7, 0.798 0.664 0.552 0.451 
9 27 1 + 0.00334 r,, 0.922 0.856 0.787 0.712 
227) 81 1 + 0.00112 r,, 0.973 0.947 0.918 0.882 
81 >81 1 1 1 1 1 
Cuadro 2. Factores para afectar la transmisividad de Darcy a los elementos del modelo MODEL 2, 
según rango de distancias y radios de Darcy. 
Valores de (1/0) para distintos radios de Darcy 
r, (m) r, (m) MODEL 2 (0) 
r,=25 m r,= 80 m 
0.25 0.8 1 +0.1183 r,, 0.2526 0.0955 
0.8 12 1 +0.05144 r,, 0.4374 0.1955 
12 5 1 +0.02736 r,, 0.5942 0.314 
3 5 1 +0.01306 r,, 0.7538 0.4889 
5 12 1 + 0.00667 r,, 0.857 0.652 
12 16 1 + 0.003625 r,, 0.917 0.775 
16 35 1 +0.00217 r,, 0.9485 0.8521 
35 40 1 + 0.001339 r,, 0.9676 0.9033 
40 100 1 +0.0008194 r,, 0.98 0.9385 
100 > 100 1 1 1 


El cuadro 2 es similar al 1, pero referido en 
este caso al modelo MODEL 2. 

La simulación se realizó para un periodo de 
24 horas con pasos de tiempo de una hora. 

Las corridas mostraron que tanto en 
régimen lineal (figura 4), como en régimen 
no lineal (figura 5), los abatimientos teóricos 
son superiores a los simulados, siendo estos 
últimos del orden del 90% de los teóricos en 
la discretización de MODEL 1 y ligeramente 
superior con la discretización de MODEL 2, 
pero del orden del 50% en el modelo regional 
(MODER 1). Esto se debe al ERD. Esta tendencia 
se mantiene a lo largo del tiempo, aunque en 
los primeros pasos de tiempo se manifiesta 
inestabilidad numérica. 

En la modelación regional de acuíferos 
siempre se ha afirmado que las cargas obteni- 


das por un proceso de simulación sólo pueden 
hacerse corresponder con los valores de las 
cargas obtenidas en los pozos de observación 
alejados de la influencia directa de estos, ya 
que los problemas locales no pueden ser 
adecuadamente reproducidos precisamente 
por el ERD. Lo anterior se evidencia en las 
figuras 4 y 5, donde se observa que a medida 
que se refina la discretización, la solución 
numérica se aproxima a la analítica. 

En estas figuras, (RL) y (RNL) significan 
régimen lineal y régimen no lineal de flujo, 
respectivamente. 

En el cuadro 3 se puede observar que 
para la distancia correspondiente al pozo de 
bombeo en el modelo regional (MODER 1), 
la diferencia entre los abatimientos obtenidos 
por la expresión analítica o teórica, y los 
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A MS A ALA - - - -| Distancia desde el 


pared 


Abatimiento RL (m) 


—4— MODER 1, RL 
—0— MODEL1,RL 


centro del pozo a la —4— MODEL 2, RL 


—"— Teórico 


Figura 4. Abatimientos teóricos y simulados en régimen lineal para un caudal de 120 1/s. 


Distancia desde el 
centro del pozo a la 
pared 


Abatimientos (m) 


-|—o MODEL 2,rp=25m,RNL |---- 


—A— Teórico, RL 


*— MODEL 1, rp=25m, RNL 


Y Teórico, rp =25 m, RNL 


—T— MODEL?2, rp= 80m, RNL 


%— MODEL 1, rp=80 m, RNL 


5 Teórico, y =80m, RNL 


Figura 5. Abatimientos simulados en régimen no lineal y su correspondencia con los teóricos, Q =1201/s. 


diferentes modelos es muy significativa y se 
van reduciendo los errores a medida que se 
refina la malla con los modelos locales MODEL 
1 y MODEL 2, siendo el error prácticamente 
independiente de la variación del caudal de 
bombeo. 


En la figura 5 se muestran los resultados 
de las corridas de MODEL 1 y MODEL 2 con 
régimen no lineal y lineal, así como la solución 
analítica en régimen no lineal. Se puede 
apreciar que el procedimiento de simular 
el régimen no lineal mediante un medio 
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Cuadro 3. Abatimientos en el pozo de bombeo en régimen lineal para los diferentes modelos y diferentes caudales. 


Abatimientos en el pozo de bombeo bajo régimen lineal en (m) 


Caudal (1/s) | Teórico a % del teórico E 1 % del teórico E 2 % del teórico 
120 225 1115 SA 2 88.89 2.1 93.33 
60 1.125 0.57 50.67 1.02 90.67 1.05 93.33 
180 Si) 17 50.96 3.06 90.67 SS 0839! 


heterogéneo equivalente es factible, y que los 
ERD se manifiestan tanto en el régimen lineal 
como en el no lineal en la misma proporción, 
en dependencia del refinamiento de la malla. 
Sin embargo, en el modelo local MODEL 1 
hay un comportamiento más errático con 
las distancias al pozo de bombeo que las 
obtenidas del modelo local MODEL 2, aunque 
los errores cometidos desde el punto de vista 
práctico indican que en la modelación regional 
una discretización del tipo de la utilizada 
en MODEL 1 para considerar los problemas 
locales en los alrededores de los pozos de 
bombeo es aceptable. 

De lo anterior se puede concluir que el flujo 
no lineal alrededor del pozo de bombeo puede 
ser perfectamente simulado como un medio 
heterogéneo mediante la reducción de la 
conductividad hidráulica de Darcy de acuerdo 
con leyes conocidas. 

De la figura 5 también se puede observar 
que los efectos de la no linealidad del flujo 
tienen su mayor incidencia en distancias muy 
próximas al pozo de bombeo (menores de 
tres metros), siendo poco significativas para 
distancias superiores a los diez metros, en 
dependencia de las propiedades del acuífero 
manifestadas en el radio de Darcy. 


Simulación de la penetración parcial de los pozos 
de bombeo 


Los efectos de la penetración parcial de 


los pozos de bombeo también pueden 
ser simulados por una variación en la 
transmisividad hidráulica, como consecuen- 
cia de la variación del espesor del acuífero. 


La variación de este espesor puede obtenerse 


partiendo de un modelo teórico que suponga 
que la línea de corriente inferior que llega 
al pozo de bombeo sigue una trayectoria 
parabólica con vértice en el radio del pozo de 
bombeo e intercepto con el impermeable a una 
distancia radial aproximadamente del orden 
del 80% del espesor del acuífero (Hernández, 
1984). Como el pozo es de pequeño diámetro, 
a los efectos prácticos se consideran las 
distancias medidas desde el centro del pozo 
de bombeo. 

De acuerdo con lo anterior, la ecuación de 
la línea de corriente sigue una parábola que 
puede ser descrita con la ecuación aproximada 
(7), brindada en (Hernández, 1984): 


Jr om 


Y=(m=h) => (0-E,) yr (7) 


Donde h, == la penetración del pozo de 
bombeo y F,==5 
(ver figura 6). 


Los espesores saturados (Esp) para las 


es el factor de penetración 


distancias utilizadas en el modelo numérico 
MODEL 1 serían: 


1=1% > Espy =h, 


m 

r=1m> Esp; a 
m 

r =3m> Esp, = Ita + (E) AS 
m 

r=9m> Esp3 =h, te 00b) 


m 
r=27m> Esp, =h, +5 (5) 43 


r=81m> Esp5 =h, +m(1-F,) 


Tecnología y 


Ciencias del Agua, vol. Il, núm. 1, enero-marzo de 2011 


Ciencias del Agua, vol. Il, núm. 1, enero-marzo de 2011 


/ 


(vy80]0u92], 


Cabrera-Estupiñán y Hernández-Valdés, Modelación del agua subterránea a escala regional con refinamiento local de la malla... 


Nivel estático 


Sp 


un 


*— hs 


Acuífero 


o | Y confinado 


Figura 6. Pozo de penetración parcial en un acuífero confinado. 


Cuadro 4. Factores de penetración parcial en función de los rangos de distancia. 


O as Factor de penetración parcial | Factor de penetración parcial | Factor de penetración parcial 
8 F (Para F =0.5) F (Para F =0.75) 
Pp Pp p Pp p 

r,>1m 0.9444 F, + 0.0556 0.528 0.7639 
lm>3m 0.8483 F, + 0.1517 0.5758 0.7879 
3m>9m 0.7372 F, + 0.2628 0.6314 0.8157 
9m>27m 0.545 FE + 0.455 0.7275 0.8638 
27m>8lm 0.2117F, + 0.7883 0.8942 0.9471 


Para cada rango de distancia se determina 
el espesor medio y éste, dividido entre el 
espesor del acuífero, determinará el factor de 
penetración parcial (F,), por el que hay que 
afectar la transmisividad hidráulica regional 
para obtener la equivalente correspondiente 
a los elementos entre dichas distancias. Así se 
tendrá lo que se observa en el cuadro 4. 

Los resultados de la simulación aparecen 
en la figura 7, donde se puede apreciar que 
los abatimientos se comportan ligeramente 
diferentes al caso anterior, ya que se reflejan 
hasta mayores distancias del pozo de bombeo 
y no como un cono tan pronunciado alrededor 
de éste. 

El abatimiento en el pozo de bombeo 
calculado por la fórmula de De Glee, citada 


por Pérez (2001), daría 2.8 y 3.89 m para 75 
y 50% de penetración, lo que confirma que 
los descensos simulados por el modelo local 
siguen siendo inferiores a los teóricos y 
que el algoritmo numérico empleado para 
representar el efecto de la penetración parcial 
es adecuado. 


Caso real de estudio “Modelo de la cuenca de 
Ariguanabo” 


El acuífero de Ariguanabo ha sido uno de los 
modelados en Cuba por su gran importancia 
para la economía del país, ya que en esta 
área se encuentran ubicados el Aeropuerto 
Internacional “José Martí” y la antigua Textilera 
“Ariguanabo”. La Cuenca Hidrográfica de 
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Modelación de la penetración parcial con MODEL 1 


Abatimiento RL 
con PP (m) 


100 


—4— Teórico en RL 


—<—  Penet. parcial 75% 
—4—  Penet. parcial 50% 


Figura 7. Variación de los abatimientos, producto de la penetración parcial del pozo de bombeo a diferentes 


distancias de éste en régimen lineal. 


Ariguanabo es una de las ocho cuencas de 
máxima prioridad del país sobre la base de su 
complejidad económica, social, ambiental, el 
grado de afectación a sus recursos naturales y 
sus características generales. En los periodos 
de intensas lluvias, este territorio se inunda 
y dificulta el funcionamiento de estas impor- 
tantes obras (Hernández et al., 2001). 

El área que se tiene en cuenta para la mo- 
delación abarca toda la cuenca Ariguanabo 
y la subcuenca Aeropuerto de la cuenca 
Almendares, en la provincia La Habana, Cuba, 
con una longitud este-oeste aproximada de 
33.5 km y unos 10 km en la dirección norte- 
sur. 

En la figura 8 se puede observar la ubica- 
ción que tiene esta cuenca en el territorio 
cubano y además se muestra la discretización 
realizada de un modelo 
AQUIMPE. 

Esta discretización regional del acuífero 


regional para 


toma la esencia de la realizada en Dilla (1988), 
donde además se llegó a calibrar el modelo, 
entiéndase obtener las propiedades K', y E por 


elementos. La diferencia fundamental radica 
en la inclusión de nuevos nodos (y por ende 
triángulos) para simular los pozos de bombeo 
CHA2, CH3 y CHA4 (ver ampliación de la figura 
8). 

Este último aspecto será importante, 
ya que en este caso de estudio se pretende 
evaluar la vinculación de la explotación real 
en cada pozo con el modelo regional. A tal 
efecto y usando la herramienta AQTRIGEO 
(ver Cabrera y Escartín, 2008), se creó un 
nuevo modelo con una discretización fina en 
cada uno de los tres pozos al estilo de MODEL 
1, quedando la zona de los pozos como se 
muestra en la figura 9, de tal forma que en este 
trabajo se emplearán dos modelos diferentes 
en cuanto a la discretización. El primero se 
llamará “Modelo Regional” y es el mostrado 
en la figura 8; el segundo es el mostrado en 
la figura 9 y se llamará “Modelo Local o de 
Explotación”. 

Para la simulación de la explotación del 
acuífero con ambos modelos, Regional y 
de Explotación, se concibe un conjunto de 
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Figura 8. Ubicación de la zona de estudio y discretización del acuífero con la numeración de nodos 


y triángulos del modelo regional. 


simplificaciones sintéticas de las caracterís- 
ticas físicas reales del sistema natural en 
estudio con fines de aplicación práctica, dentro 
de las que está la de acuífero confinado con 
simulación bidimensional en planta, régimen 
impermanente con medio continuo poroso 
isotrópico y heterogéneo. 

La simulación se realizará en un año divi- 
dido en doce intervalos de tiempos de treinta 
días cada uno. Los datos referentes al modelo 
fueron tomados de Dilla (1988) y procesados 
para adaptarlos al presente trabajo. 

La recarga, la explotación al acuífero y 
los niveles referidos al estado inicial para el 
periodo en estudio también fueron tomadas 
de Dilla (1988). En este sentido, es importante 
destacar que se realizó una redistribución de 
las extracciones para garantizar que de los 
nodos 75, 79 y 82, asignados a los pozos CH2, 


CH4 y CH3 del modelo regional, se extrajeran 
caudales de 200 1/s en los casos de los dos 
primeros pozos, y 220 1/s en el caso del CH3. 
Estos son los valores reales posibles a extraer 
en dichos pozos. De igual forma, estos caudales 
fueron asignados a los nodos correspondien- 
tes del Modelo de Explotación, con el fin de 
realizar comparaciones entre ambos modelos. 

Además de la simulación numérica con 
AQUIMPE, se decidió realizar una evaluación 
de las expresiones teóricas en cada pozo tanto 
de flujo lineal como no lineal y aplicar el 
procedimiento propuesto para introducir el 
efecto de la no linealidad mediante la variación 
gradual de la K,, de los elementos vecinos a los 
pozos. 

Los resultados pueden ser observados 
en las ilustraciones 10, 11, 12, las cuales son 
comentadas a continuación. 
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Figura 9. Discretización del “Modelo de Explotación”, donde se ubican mallados finos en los pozos CH2, CH3 y CH4. 


Abatimientos alrededor del pozo CH3 obtenidos por varios modelos 


Niveles (m) 


Modelo explotación, RL 
Teórico, RL 
Teórico, RNL 


Regional, RL 
Modelo explotación, RNL 


Distancia (m) 


Figura 10. Conos de abatimientos en la vecindad del pozo CH3 para el tiempo + 12, considerando varios 
modelos para su determinación. 


En la figura 10 se muestra un conjunto de La serie “Regional, RL” está definida con 
series obtenidas para el pozo CH3 en el tiempo un solo valor de cota del agua en el pozo. 
12 o final del último mes de explotación, para Este valor es 38.09 m y fue obtenido mediante 
una distancia de hasta 10 m desde el eje del la simulación con AQUIMPE, empleando el 
pozo de explotación. modelo Regional mostrado en la figura 8, 
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donde el pozo de referencia es modelado con 
un solo nodo y a su alrededor no existe un 
mallado fino. 

La serie RL” 
evaluando la parte lineal o darciana de la 
expresión (1), utilizando una T,,= 2 500 m?/ d, 
un E=0.15, y un caudal de 220 1/s durante un 
tiempo de cuatro meses de bombeo continuo, 


“Teórico, fue obtenida 


tiempo suficiente para lograr la estabilidad. 
Además, esta expresión fue evaluada para 
diferentes radios medidos desde el eje del 
pozo. 

La serie “Modelo explotación, RL” fue 
obtenida aplicando AQUIMPE, pero en este 
caso se utilizó la discretización del modelo 
de explotación, en donde todos los triángulos 
que se encontraban dentro del elemento 
definido por los nodos 68, 232 y 578 del 
modelo de explotación tienen las mismas 
propiedades K,, m y E, que el elemento 
definido por los nodos 68, 76 y 90 del modelo 
regional. 

Hasta el 
diferencia de abatimientos de alrededor de 


momento cabe destacar la 
siete metros entre la serie “Regional, RL” con 
las obtenidas en “Teórico, RL” y en “Modelo 
explotación, RL”, lo cual demuestra que para 
tratar de simular con menor error los niveles 
piezométricos en los pozos de bombeo no basta 
con considerar un esquema o discretización 
regional, sino que se debe realizar un mallado 
fino alrededor del pozo para evitar en la mayor 
medida posible los Errores de Refinamiento de 
la Discretización (ERD). 

Por si esto fuera poco, cuando se evalúa 
la expresión (1) en su totalidad, teniendo en 
cuenta un radio de Darcy de 80 m, se obtiene 
la serie “Teórico, RNL”, en la cual se puede 
observar que en los alrededores del pozo 
ocurren descensos importantes de los niveles 
debidos a la no linealidad del flujo en esa zona. 
Por esta razón hay un aumento considerable 
de las pérdidas energéticas, lo cual provoca 
que bajo estas condiciones la diferencia entre 
niveles en el pozo obtenidos por los modelos 
“Regional, RL” y “Teórico, RNL” sea de 
aproximadamente unos 23 m. 


Para obtener la variación gradual de la K,, 
de los elementos vecinos a los pozos y así tener 
en cuenta el efecto de la no linealidad del flujo 
en el modelo numérico, se empleó la corrección 
brindada en la expresión (6). Se tomaron los 
valores necesarios del cuadro 1, aunque se 
debieron realizar algunos ajustes, ya que por 
ejemplo, estos tres pozos modelados tienen un 
radio de 30 cm. 

Una vez que se afectaron las propiedades 
de los triángulos vecinos a cada pozo se realizó 
una simulación con AQUIMPE, teniendo 
en cuenta estos cambios. En las figuras 10 y 
11 se puede observar que existe una buena 
coincidencia entre las series “Teórico, RNL” y 
“Modelo Explotación, RNL”. 

Estos resultados tienen un valor práctico 
muy importante: se tiene un modelo numérico 
que considera no sólo los efectos regionales 
del acuífero, sino que como producto de 
esta recalibración con sólo los elementos 
más próximos a los pozos de explotación se 
pueden obtener los niveles más próximos a 
los reales en esta zona bajo cualquier régimen 
de recarga-explotación. Esto permite evaluar 
con mayor certeza posibles fenómenos 
locales, como la intrusión salina, problemas 
de contaminación de los pozos, etcétera. 

En la figura 11 se observan los mismos 
resultados presentados en la figura 10, pero 
se aumenta la distancia a la que se evalúan 
los modelos, llegando hasta una longitud de 
81 m, medidos desde el eje del pozo CH3. En 
estas dos gráficas, así como en la mostrada 
por la figura 4, se puede observar y corroborar 
que el efecto de la no linealidad se desarrolla 
en distancias cercanas a los pozos de bombeo, 
unos 10 m, aproximadamente. A partir de 
esta distancia, todos los modelos comienzan 
a coincidir. 

Finalmente se muestra la figura 12, que 
contiene la planta a), donde se muestran los 
niveles del agua subterránea para el tiempo 
12 en la zona donde se encuentran los pozos 
analizados. Este mapa fue obtenido mediante 
la simulación, teniendo en cuenta el efecto de 
no linealidad, aunque éste no puede apreciarse 
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Abatimientos versus longitudes medidas desde el radio del pozo CH3 


—— Teórico, RL 
aus Pared delpozo | ___| —+— Modelo explotación, RL 


—— Teórico, RNL 


—-— Modelo explotación, RNL 


Abatimiento (m) 


Figura 11. Abatimientos en el pozo CH3 obtenidos por varios modelos y evaluados desde el radio del pozo 
hasta una distancia de 81 m. 


347 500 


347 000 


Pozo CH4 
Pozo CH2 Pozo CH3 


b) 


36 a) Mapa de niveles del agua subterránea correspondiente 
al tiempo núm. 12 y concentrado en la zona de los pozos 
34 CH2, CH3 y CH4. 
b) Vista tridimensional de los conos de abatimiento 
producidos en la vecindad de los pozos. 


345500 346000 346500 347 000 347500 348 000 348 500 349 000 a) 


Figura 12. Mapa donde se muestran los conos de abatimiento en los pozos de bombeo CH2, CH3 y CH4. 


en toda su magnitud, ya que la interpolación 14 m en el pozo CH3 no se pueden apreciar. 
fue realizada con celdas regulares de 10 m Simplemente es un problema de escala. 
por 10 m. Por tal razón se muestran las zonas También aparece la imagen b), donde se 
de depresión en color azul. Sin embargo, los muestran de forma tridimensional los conos 
detalles de niveles que llegan a valores de de abatimiento en cada uno de los pozos y con 
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una vista inferior, ilustrando las posibilidades 
gráficas que presenta esta versión del modelo 
de flujo AQUIMPE desarrollado sobre Matlab. 


Conclusiones 


Con este ensayo numérico se ha demostrado 
que un refinamiento de la discretización, 
como el utilizado en el modelo MODEL 1 
alrededor de cada pozo de bombeo y su acople 
al modelo regional, es un adecuado algoritmo 
para simular los efectos locales que provocan 
abatimientos significativos en los pozos de 
bombeo y sus proximidades, los cuales no 
se manifiestan en el modelo regional, ya sea 
como: 


e Producto de la presencia de un flujo lineal 
en todo el acuífero y en los alrededores del 
pozo de bombeo. 

e Efectos de la no linealidad del flujo emplean- 
do un medio heterogéneo equivalente, 
afectando la transmisividad calibrada en el 
modelo regional por factores de acuerdo con 
rangos de distancias del pozo de bombeo. 

e Efecto de la penetración parcial de los pozos 
de bombeo. 


Se brinda una estrategia O procedimiento 
sencillo y eficaz para tener en cuenta efectos 
locales en los pozos de explotación, como 
pueden ser la presencia de flujo no lineal 
y la penetración parcial del pozo en el 
acuífero. Este procedimiento está basado en 
la recalibración gradual de las propiedades 
de los elementos vecinos al pozo de bombeo, 
teniendo en cuenta un mallado fino alrededor 
de cada pozo de explotación analizado. 
Estos elementos enriquecen las posibilidades 
que brinda la simulación matemática para 
el mejor entendimiento de fenómenos tan 
complejos como el movimiento de las aguas 
subterráneas. 

Se aplican los procedimientos propuestos 
en el caso de estudio del acuífero de Arigua- 
nabo, utilizándose tres pozos de bombeo, 
en los cuales se concentró la explotación. 


Se aplicaron los Modelos Regional y de 
Explotación para analizar el comportamiento 
de los niveles tanto en régimen lineal como no 
lineal. 
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Abstract 


CABRERA-ESTUPIÑÁN, E. € HERNÁNDEZ-VALDÉS, A. Regional groundwater 
modeling with a local mesh refinement. Algorithm presentation and validation. Water 
Technology and Sciences, formerly Hydraulic engineering in Mexico (in Spanish). Vol. 
II, No. 1, January-March, 2011, pp. 65-82. 


Groundwater flow simulation at a regional scale cannot accurately reproduce significant 
drawdowns in the well fields, because in general these models have been used to assess the 
hydrodynamic behavior of the aquifer system and hydraulic system effects at such scale. 
Presently, work is being done on the validation of an algorithm which, after starting with a 
regional model, is able to reproduce the local behavior of levels in wells located near selected 
pumping wells. The pumping wells are those which require a systematic control as an 
operating object in real time by means of automatic control techniques. By means of validation 
with analytic solutions, the importance of the refinement of the discretization is demonstrated, 
even when using the quadratic triangle in the Finite Element Method which is the one used in 
the AQUIMPE technology, in its recent Windows version. The proposed algorithm allows to 
reproduce the effects of the non-linearity of flow in the vicinities of the pumping wells and the 
effects of partial penetration by introducing suitable modifications to aquifer transmissivity in 
the elements near to the pumping well according to the distance ranges used by the numeric 
algorithm. 


Keywords: mathematical modelling, aquifers, water resources, geographic information 
systems, finite element method, mesh refinement. 
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